library(haven)
library(ggplot2)
library(tidyr)
library(dplyr)
library(sqldf)
library(readxl)

setwd("C:/Users/jonves/Dropbox/Papers/SSP-Instability-Opinion/replication_data/")

ssp <- readRDS("data/ssp_long.rds")

#### Forecast data ####

models <- c("OECD Env-Growth", "IIASA GDP", "IIASA-WiC POP")
scenarios <- c("SSP1", "SSP2", "SSP3", "SSP4", "SSP5")
years <- c(2010, 2100)

variables <- c("region", "year", "scenario", "model", "value")

df <- ssp[which(ssp$model %in% models 
                 & ssp$scenario %in% scenarios
                 & ssp$year %in% years), variables]

df <- spread(df, model, value)


df$OECD <- df$`OECD Env-Growth`/df$`IIASA-WiC POP`*1000
df$IIASA <- df$`IIASA GDP`/df$`IIASA-WiC POP`*1000
df <- df[c("region", "year", "scenario", "OECD", "IIASA")]

df <- gather(df, model, gdpcap, -region, -year, -scenario)
df$model <- paste(df$model, df$scenario)
df <- df[c("region", "year", "model", "gdpcap")]
df <- spread(df, year, gdpcap)

names(df) <- c("region", "model", "gdpcap_2010", "gdpcap_2100")

df$avg_growth20102100 <- log(df$gdpcap_2100/df$gdpcap_2010)
df$x <- log(df$gdpcap_2010)

#### Historical data ####
# SSP historical data useless for doing growth disruptions as it only is reported at 5 year intervals.
wdi <- filter(ssp, scenario %in% c("History - World Bank (WDI)", "History - WPP2010"))



pwt <- read_dta("data/pwt90.dta")
pwt <- pwt[which(!is.na(pwt$rgdpna)),]

pwt$gdpcap <- pwt$rgdpna / pwt$pop
growth <- function(x){x/lag(x)-1}

pwt <- group_by(pwt, countrycode) %>%
  arrange(year) %>%
  mutate(grwt = growth(gdpcap)) %>%
  summarize(ngrwtr = mean((grwt < 0), na.rm = T))


# Merge with SSP projections
df <- left_join(df, pwt, by=c("region" = "countrycode"))

df$x <- df$ngrwtr
df$y <- df$avg_growth20102100


lm_eqn <- function(y, x){
  m <- lm(y ~ x)
  eq <- substitute(italic(y) == a + b * italic(x), 
                   list(a = as.numeric(round(coef(m)[1], digits = 1)), 
                        b = as.numeric(round(coef(m)[2], digits = 1)), 
                        r2 = as.numeric(round(summary(m)$r.squared, digits = 2))))
  as.character(as.expression(eq))
}

getr2 <- function(y, x){
  m <- lm(y ~ x)
  r2 = format(summary(m)$r.squared, digits = 3)
}

df <- group_by(df, model) %>%
  mutate(label=lm_eqn(y, x))


p1 <- ggplot(df, aes(x=x, y=(y / (2100-2010)), label=label)) + geom_point(size=0.5) +
  geom_smooth(method = "lm", se=FALSE, color="red", formula = y ~ x) +
  facet_wrap(~model, nrow=2, ncol=5) +
  geom_text(x=0.35, y=6, size=4, parse=TRUE) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  theme_bw() + xlab("Share of years with negative growth (PWT 9.0, 1950-2014)") + ylab("Average yearly growth rates (SSP, 2010-2100)") +
  theme(strip.text.x = element_text(size=14),
        axis.title = element_text(size=14),
        axis.text = element_text(size=14))

ggsave("results/figure3.tiff", plot=p1, width=7.5, height=6, dpi=300, compression="zip")


